Modelling Deformation of Lake Sediments from Coring and Extrusion
Dewey Dunnington
July 9, 2016
library(ggplot2)
library(reshape2)
library(dplyr)Introduction
Lake sediments get extruded and deformed like this:
plot(c(-1, 1), c(0, 10), ylim=c(10, 0), pch=".", asp = 1)
rect(-1, 10, 1, 0)
xvals <- seq(-1, 1, 0.02)
for(d0 in seq(0, 9, 0.2)) {
lines(xvals, sapply(xvals, function(x) x^2+d0))
}
rm(d0, xvals)If we take a slice through this (as we would while extruding), what is the distribution of original depths? Does the slice size or degree of deformation affect this original depth distribution? Is it possible to correct for deformation to “sharpen” the data we infer from cores?
plot(c(-1, 1), c(0, 10), ylim=c(10, 0), pch=".", asp = 1)
rect(-1, 10, 1, 0)
xvals <- seq(-1, 1, 0.02)
for(d0 in seq(0, 9, 0.2)) {
lines(xvals, sapply(xvals, function(x) x^2+d0))
}
rect(-1, 2, 1, 3, lwd=2)
rm(d0, xvals)Creating a Model
Imagine we have a 2 cm wide core where we can model deformation as d = d0 + r^2 (where r is the distance to the centre of the core barrel and d0 is the original depth prior to deformation). It follows that we can then model the original depth of the sediment at a position (d, r) as d0(d, r) = d - r^2. Rasterizing this, we can make a model of actual and original depths given a function d0(d, r)
deformation_model3d <- function(slicesize, d0, width=2, cellsize=0.025) {
coords <- expand.grid(x=seq(-width/2, width/2, cellsize),
y=seq(-width/2, width/2, cellsize),
z=seq(-slicesize/2, slicesize/2, cellsize))
coords$r <- sqrt(coords$x^2 + coords$y^2)
coords <- coords[coords$r <= width/2,]
coords$d0 <- d0(coords$z, coords$r)
return(coords)
}This function takes a slicesize (for example, 1 cm) and a function d0 that returns a the original depth prior to deformation given a depth (d) and a radius from the centre of the core (r). The cellsize parameter specifies the resolution at which to create the model (low values for higher resolution) and the width specifies the interior diameter of the core barrel. The d values in the model are relative to the centre of the slice. For a single function d0 and slice size, we can then plot a cross-section and histogram of the depths prior to deformation.
model <- deformation_model3d(1, function(d, r) d - r^2)
# overhead view (same for all sizes of slice)
ggplot(model %>% filter(z==0), aes(x=x, y=y)) +
geom_raster(aes(fill=d0)) +
scale_fill_gradient2() + scale_y_reverse() +
coord_fixed()
# side on view
ggplot(model %>% filter(y==0), aes(x=x, y=z)) + geom_raster(aes(fill=d0)) +
scale_fill_gradient2() + scale_y_reverse() +
coord_fixed()
# density histogram
ggplot(model, aes(x=d0)) + geom_density() 
This histogram seems to indicate that the average actual depth (NA) is slightly smaller than the centre of the slice.
Effect of degree of deformation and slice size
We can also vary these parameters to model the effect of slice size and degree of deformation.
# dplyr doesn't like passing functions or evals so we are going to wrap this
# so we can pass a single number
deformation_model_bydegree <- function(slicesize, degree, cellsize=0.025, width=2) {
d0 <- function(d, r) d - abs(r^degree)
deformation_model3d(slicesize, d0, width, cellsize)
}
# create the models in long form so we can plot with facets
all <- expand.grid(slicesize=seq(0.2, 1, 0.2), degree=2:6) %>%
group_by(slicesize, degree) %>%
do(deformation_model_bydegree(.$slicesize, .$degree))
# overhead view (same for all sizes of slice)
ggplot(all %>% filter(z==0, slicesize==all$slicesize[1]), aes(x=x, y=y)) +
geom_raster(aes(fill=d0)) +
scale_fill_gradient2() + scale_y_reverse() +
coord_fixed() + facet_grid(slicesize ~ degree)
# side on view
ggplot(all %>% filter(y==0), aes(x=x, y=z)) + geom_raster(aes(fill=d0)) +
scale_fill_gradient2() + scale_y_reverse() +
coord_fixed() + facet_grid(slicesize ~ degree)
# density histograms
ggplot(all, aes(x=d0)) + geom_density() + facet_grid(slicesize ~ degree)
# we can also create a table summarizing some things
all_summary <- all %>% summarise(mean_depth=mean(d0))knitr::kable(dcast(all_summary, degree~slicesize), digits=2)## Using mean_depth as value column: use value.var to override.| degree | 0.2 | 0.4 | 0.6 | 0.8 | 1 |
|---|---|---|---|---|---|
| 2 | -0.50 | -0.50 | -0.50 | -0.50 | -0.50 |
| 3 | -0.40 | -0.40 | -0.40 | -0.40 | -0.40 |
| 4 | -0.33 | -0.33 | -0.33 | -0.33 | -0.33 |
| 5 | -0.29 | -0.29 | -0.29 | -0.29 | -0.29 |
| 6 | -0.25 | -0.25 | -0.25 | -0.25 | -0.25 |
Looks like for the even powers this is pretty easy. The weighted average depth of a given slice is always -1/(degree/2). This means the average depth is always slightly less than (25 to 50%) the value stated, regardless of slice thickness.
rm(all, all_summary, model, deformation_model_bydegree)An assessment of actual deformation in cores
Method: ImageJ to scale pictures of deformation in cores where deformation does occur.
# load
deformations <- read.csv("deformed cores/deformed_cores_summary.csv")
# calc the x0 and y0 to use as the base coordinates (to transform x and y to r and d)
layers <- deformations %>%
group_by(layercode) %>%
summarise(x0=mean(range(x)), y0=min(y))
deformations <- merge(deformations, layers, by="layercode")
deformations$r <- deformations$x-deformations$x0
deformations$d <- deformations$y-deformations$y0
deformations <- deformations %>% select(layercode, core, layer, r, d)
# plot cores with quadratic smoothing
ggplot(deformations, aes(x=r, y=d)) +
geom_point(aes(col=factor(layer))) +
stat_smooth(method=lm, formula=y ~ poly(x, 2, raw=TRUE), se=FALSE) +
scale_y_reverse() + facet_wrap(~core, scales="free")
# calculate quadratic models
modelparams <- deformations %>%
group_by(core) %>%
do(data.frame(t(data.frame(lm(data=., formula=d ~ poly(r, 2, raw=TRUE))$coefficients))))
names(modelparams) <- c("core", "c", "b", "a") # in y=ax^2 + by + c
modelparams <- modelparams %>% select(core, a, b, c)knitr::kable(modelparams, digits=2)| core | a | b | c |
|---|---|---|---|
| longlake_pc1 | 0.51 | 0.10 | 0.22 |
| menounos_cheak1 | 0.18 | 0.02 | 0.01 |
| menounos_cheak2 | 0.33 | 0.11 | 0.00 |
| sl2 | 0.11 | 0.07 | -0.01 |
| suzie1 | 0.32 | 0.16 | -0.18 |
| whistler_gc4 | 0.16 | -0.02 | -0.02 |
| whistler_gc8 | 0.10 | 0.09 | 0.00 |
A more realistic model
It appears a quadratic is a good fit for experimental deformation data (where cores are deformed, which not all are). Coefficients for x^2 range from 0.1 to 0.6. Core barrel diameters range from 2 to 5 inches (5 to 13 cm), so we will use these end members as our test cases.
Effect of slice size
deformation_model_wrapper <- function(slicesize, coeff, width) {
deformation_model3d(slicesize, d0=function(d, r) d - coeff * r^2, width = width, cellsize = width/100)
}
# create the models in long form so we can plot with facets
all <- expand.grid(slicesize=seq(0.2, 1, 0.2), coeff=c(0.1, 0.3, 0.5)) %>%
group_by(slicesize, coeff) %>%
do(deformation_model_wrapper(.$slicesize, .$coeff, 7))
# overhead view (same for all sizes of slice)
ggplot(all %>% filter(abs(z)==min(abs(z)), slicesize==all$slicesize[1]), aes(x=x, y=y)) +
geom_raster(aes(fill=d0)) +
scale_fill_gradient2() + scale_y_reverse() +
coord_fixed() + facet_grid(slicesize ~ coeff)
# side on view
ggplot(all %>% filter(abs(y)==min(abs(y))), aes(x=x, y=z)) + geom_raster(aes(fill=d0)) +
scale_fill_gradient2() + scale_y_reverse() +
coord_fixed() + facet_grid(slicesize ~ coeff)
# density histograms
ggplot(all, aes(x=d0)) + geom_density() + facet_grid(slicesize ~ coeff)
Effect of barrel width
# create the models in long form so we can plot with facets
all <- expand.grid(width=c(5, 7, 13), coeff=c(0.1, 0.3, 0.5)) %>%
group_by(width, coeff) %>%
do(deformation_model_wrapper(1, .$coeff, .$width))
# side on view
ggplot(all %>% filter(y==0), aes(x=x, y=z)) + geom_raster(aes(fill=d0)) +
scale_fill_gradient2() + scale_y_reverse() +
coord_fixed() + facet_grid(width ~ coeff)
# density histograms
ggplot(all, aes(x=d0)) + geom_density() + facet_grid(width ~ coeff)
Smoothing effects on stratigraphic data
The distribution (density) acts as a smoothing filter on geochem data.
Conclusions
There is a limit to how small extrusion intervals can get based on deformation (flat topped distributions are bad because they don’t have a ‘mode’ depth that they represent!).